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George  D.  Ashton 

INTRODUCTION 

The  use  of  air  bubbler  systems  to  suppress  ice 
formation  by  inducing  a flow  of  water  against  the 
underside  of  an  ice  cover  is  a commonly  used  tech- 
nique. The  analysis  of  line  source  bubblers  for  such 
purposes  has  recently  (Ashton  1974, 1978)  progressed 
to  the  extent  that  it  is  possible  to  simulate  and  pre- 
dict the  performance  of  these  bubblers,  and  that  anal- 
ysis has  been  validated  against  field  and  laboratory 
data  (Ashton  197':,  1978).  This  report  presents  a 
parallel  analysis  for  point  source  bubbler  systems,  i.e., 
the  suppression  of  ice  that  results  from  a single-point 
discharge  of  air  bubbles.  The  installations  of  such 
systems  would  be  useful,  for  example,  in  the  protec- 
tion of  individual  piles  of  a multipilc  dock  installation. 


OUTLINE  OF  ANALYSIS 

This  study  begins  with  an  analysis  of  the  plume 
induced  by  the  rising  stream  of  bubbles  emanating 
from  an  orifice  submerged  in  the  water  and  uses  the 
results  of  Kobus  (1%8).  It  next  determines  the  heat 
transfer  coefficient  using  the  plume  parameters  at  the 
point  of  impingement  of  the  plume  on  the  underside 
of  the  ice  cover  by  analogy  with  empirical  results  of 
Gardon  and  Akfirat  (1966)  for  impinging  axisym- 
metric  air  jets.  It  thus  determines  ice  melting  and 
thermal  depletion  by  use  of  a simplified  energy  bud- 
get calculation.  Finally,  it  applies  the  resulting 
analysis  to  the  practical  case  of  varying  winter  temper- 
atures by  a quasi-steady  stepwise  solution  utilizing 
daily  temperatures.  An  example  simulation  is  pre- 
sented. The  FORTRAN  computer  program  for  the 
simulation  is  given  in  the  Appendix. 


PLUME  ANALYSIS 


The  air  discharge  rate  Qq  (m*  s'')  from  an  orifice 
of  diameter  d (m)  is 


Oo-c,f- 


(1) 


where  is  a discharge  coefficient  (on  the  order  of 
0.6;  see  e.g..  Rouse  1946),  ^ is  the  pressure  dif- 
ference across  the  orifice,  and  p,  is  the  air  density 
inside  the  bubbler  line.  Typical  orifice  diameters 
used  in  existing  installations  are  on  the  order  of  1 mm. 
The  pressure  difference  ^ is 

Ap  = (2) 


where  is  the  pressure  inside  the  supply  line, 

p^gH  is  the  hydrostatic  pressure  at  the  submergence 
depth  H in  water  of  density  p^,  and  g is  the  gravita- 
tional constant. 

Using  the  results  of  Kobus  (1968),  the  centerline 
water  velocity  U^(x)  (m  s'*  )at  distance  x above  the 
orifice  is  given  by 

‘ C (x+Xq)  L J (3) 

where  c is  the  rate  of  linear  spread  of  the  plume,  Xq 
is  an  empirical  coordinate  correction  to  account  for 
various  near-orifice  effects,  is  the  atmospheric 
pressure,  H*  = H-^P^tmlp^g  (for  sea  level  conditions 
H*  = /y+10.3  m),  and  is  the  mean  rising  speed  of 
the  bubbles.  These  are  further  illustrated  in  Figure  1 . 
Both  and  c were  found  by  Kobus  (1%8)  to  be 
weak  functions  cf  (?q  and  a fit  of  these  data  yields 

(5) 

where  Cj  = 0.1 52  m'® ’ s*  ' ’ and  Cb  = 1 .83  m®  * * 
s'*  ■*  ’ . The  plume  analysis  of  Kobus  used  a Gaussian 
distribution  of  vertical  water  velocity  of  the  form 


^ (*/  = exp 

U,{x) 


2c^  (x+Xq)^ 


(6) 


where  r is  the  radial  coordinate  from  the  plume 
centerline.  The  total  volume  flux  (x)  is  then 
given  by 
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Figure  J.  Definition  sketch.  Figure  2.  Centerline  velocity  of  plume  as  a function  of 

orifice  discharge  and  submergence  depth  of  orifice. 


M (m) 

Figure  3.  Induced  flow  of  water  at  Impingement  as  a function 
of  air  discharge  rate  and  submergence  depth  of  orifice. 
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Figure  4.  Diameter  of  impinging  plume  as  a function  of 
air  discharge  rate  and  submergence  depth  of  orifice. 


Q*(x)  = 2nt/,  (x)cMx+Xo)V  (7)  6=(/y+Xo)C, 


Appropriate  substitution  for  c and  Uy^  into  eqs  3 
and  7 then  yields,  respectively, 

C,  (x+Xq)  [ itp^  Cy,  J (8 


(X)  = 2C,  (x+Xq)  (?o® 


r-^en»^'og.l1-(W/V*)n°-^  (9) 

L J 

Defining  (x)  and  = (x)  at  X = H 

yields 

u r-/>,unios,n-w^*)iT-^ 

C,(W+Xo)[  rrp^Cb  J(10) 

= 2C,  (W+Xo)  Oo®-*’* 

(11) 

L Cb  J 


In  Figures  2,  3,  and  4 are  shown,  respectively,  the 
values  of  U^y^,  C?„h>  ^ functions  of  air  dis- 

charge rate  C?o  and  submergence  depth  H. 


HEAT  TRANSFER  ANALYSIS 

The  heat  transfer  analysis  consists  of  determining 
the  temperature  of  the  plume  impinging  on  the 
underside  of  the  ice  cover,  estimating  the  heat  trans- 
fer coefficient  at  the  underside  of  the  cover,  perform- 
ing a simplified  analysis  of  the  thickening  (or  thinning) 
of  the  cover,  and  depleting  the  thermal  reserve  as  a 
consequence  of  the  heat  transfer  to  the  cover. 

Temperature  of  impinging  plume 

Since  the  water  bodies  in  which  bubbler  systems 
are  installed  are  seldom  isothermal,  it  is  necessary 
to  evaluate  the  entrainment  of  water  at  different 
levels  above  the  point  of  air  discharge  to  arrive  at  the 
‘mixed’  temperature  at  impingement.  That  is,  the 
impingement  temperature  referenced  to  the 
freezing  point  r„,  is  given  by 


VwH  0 


and  the  width  b at  x = //  is 


3 


r 


I 

I 


If 


1 

I 


i 

i 

I 


where  <y  (?*  {x)/dx  is  the  entrainment  rate  and  r^(jc) 
is  the  water  temperature  as  a function  of  x above  the 
air  discharge  point.  For  a vertically  uniform  ambient 
water  temperature  the  impingement  tempera- 
ture = ^wA-  temperature  profiles  may 

easily  be  integrated  using  eq  1 3 to  arrive  at 

Heat  transfer  coefficient 

Gardon  and  Akfirat  (1966)  found  the  heat  trans- 
fer coefficient  associated  with  an  axisymmetric  im- 
pinging air  jet  (in  air)  to  be  correlated  by 

(14) 


where  the  Nusselt  number  averaged  over  an  area  of 
diameter  2 Tq  defined  by 


Nu. 


= '"O 


(15) 


and  is  the  average  heat  transfer  coefficient  and  k 
is  the  thermal  conductivity  of  the  air.  The  associated 
Reynolds  number  is  defined  by 


(16) 


where  is  the  axial  velocity  of  the  air  jet  and  p and 
p are  the  density  and  viscosity  of  air. 

By  analogy  with  other  heat  transfer  results  (see, 
e.g.,  Rohsenow  and  Choi  1961 ),  we  may  convert  eq 
14  to  a more  general  relationship  by  introducing  the 
Prandtl  number  dependence  in  the  form 


Nu  <x  Pr’^^ 

where  Prandtl  number  Pr  is  defined  by 


(17) 


(18) 


The  heat  transfer  coefficient  itr  = b will  be  used 
as  a reference  value  with  which  to  normalize  the 
radial  variation  of /).  Thus 


h^,  = 2.08- 


k U. 


cH 


035  fc-0.4S 


.,035 


(20) 


where  v = p/p  and 

h(r)=h^(rlb)-°*^  (21) 

for  /■  > 6,  For  r<b,v/e  simply  fit  a parabola  to 
eq  21  such  that  dh(r)ldr  = 0 at  r = 0 and  has  the 
same  slope  and  magnitude  as  eq  21  at  r = b.  Hence, 
for  r < 6 


h(r)  = /tb  1 1 .225-0.225  (r/bf  \ . (22) 

Variation  of  /?(,  as  a function  of  submergence  depth 
and  air  discharge  is  presented  in  Figure  5.  The 
variation  oi  h(r)  .normalized  by  /?(,  is  presented  in 
Figure  6 as  a func  rion  oi  rib  for  r>  b. 

Melting  of  the  ice  cover 

The  actual  melting  of  the  ice  is  governed  by  the 
heat  balance  at  the  water/ice  interface,  given  by; 

(«) 

where  is  the  rate  of  heat  conduction  through  the 
ice,q^  = ft(7'^H‘7’m)  t^e  rate  of  heat  transfer  to 
the  undersurface  of  the  ice,Pj  is  the  mass  density  of 
the  ice,  X is  the  heat  of  fusion  of  the  ice,  and  dq  jdt 
is  the  rate  of  change  of  the  ice  thickness. 

The  model  used  for  Pj  is  one-dimensional  steady- 
state  heat  conduction,  which  assumes  a linear  vari- 
ation in  temperature  through  ice  thickness  (and  snow 
thickness)  together  with  an  estimate  of  the  transfer 
coefficient  through  the  air  boundary  layer.  Thus 


and  Cp  is  the  specific  heat  of  the  fluid.  Assuming  Pr 
= 0.7  for  air  and  Pr  = 1 3.6  for  water  at  O'C  enables  eq 
14  to  be  transformed  to  apply  to  water  flow  in  the 
form 


= 2.08  Re*®-**  (19) 

where  the  Reynolds  number  is  now  that  for  water; 
tha>  is.  Re*  = <7^^  ftg  p*/p. 


9i  = 


(24) 


where  is  the  top  surface  temperature  of  the  ice 
(or  snow,  if  present),  and  bj  and  are  the  thermal 
conductivity  of  ice  and  snow.  In  general,  Tj  is 
not  the  ambient  air  temperature.  Within  the  con- 
text of  the  steady-state  assumptions  above,  it  is 
reasonable  to  represent  the  heat  transferred  through 
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Figure  5.  Variation  of  as  a function  of  orifice  discharge  and 
submergence  depth  of  orifice. 
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Figure  6.  Variation  of  local  heat  transfer  rate  as  a function  of  radial 
distance  from  centerline  of  Impingement. 


the  air  boundary  layer  in  the  form 


- ( 

Uh 


a 


(25) 


where  l/hj  represents  a thermal  resistance  due  to  the 
air  boundary  layer.  Assuming  values  in  eqs  24  and 
25  are  equal  (equivalent  to  a scries  representation  of 
the  thermal  resistances),  then  f,  may  be  eliminated 
and 


(t?i/<fi)+(»Js/^s)+(1/^a) 


(26) 


The  difficulty  is  to  obtain  a reasonable  estimate  of 
)obson  (1973)  examined  the  energy  budget 
terms  and  found  that  a good  approximation  for  ft, 
from  a water  surface  is  of  the  form 

= 3.4+4.4  (27) 
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Figure  7.  A ir  temperature  record  used  in  simulation  example. 
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Figure  8,  Variation  of  open  water  area  for  simulation 
example. 


where  is  the  wind  velocity  (m  s'* ) and  h,  is  the 
heat  transfer  coefficient  in  W . 

In  the  following  simulation  example,  we  will 
arbitrarily  take  = 24  W m'^°C’  corresponding 
approximately  to  a windspeed  of  4.5  m s"' . Although 
is  derived  for  open  water,  we  will  assume  it  is  also 
applicable  to  the  ice/air  interface.  Similarly,  the 
heat  loss  from  open  water  above  the  bubbler  uses 


eq  26  with  = T.^  and  the  same  value  of  h^.  It 
is  recognized  that  a more  detailed  energy  budget 
approach  could  be  used,  but  the  detail  is  considered 
inappropriate  in  view  of  other  uncertainties  in  the 
analysis  and  the  lack  of  detailed  input  data  other 
than  daily  maximum  and  minimum  air  temperatures 
usually  available.  = 24  W m*^  °C"' , incidentally, 
is  equivalent  to  the  thermal  resistance  of 
approximately  0.0%  m of  solid  ice. 


( Olitooc*  from  ^ tml 


Figure  9.  Ice  thickness  profile  on  28  lanuary 
of  simulation  example. 


Simulation  example 

Figure  7 shows  a daily  average  air  temperature 
record  taken  from  a midwestern  location.  It  shows 
the  patterns  of  alternating  very  cold  periods  with 
somewhat  warmer  pericxis  which  arc  typical  of 
winter  periods.  F igure  8 shows  the  results  of  using 
that  record  in  the  present  simulation  as  a plot  of 
radius  of  open  water  extent  with  time  for  the  param- 
eters of  = 0.0005  m’  s'*  (^1  ftVmin),  = 
0.2‘’C.  //  = 5 m,  and  /r^  = 25  . 

Figure  8 shows  that  the  extent  of  open  water 
responds  to  the  air  temperature  variation  by  con- 
tracting during  cold  periixJs  and  expanding  during 
warmer  peritxis.  While  the  ice  cover  is  predicted  to 
free/e  over  the  bubbler  during  very  cold  days,  the 
thickness  is  considerably  reduced  from  that  which 
would  exist  if  the  bubbler  had  not  been  present,  as 
shown  in  the  simulated  profile  of  ice  thickness  for 
28  January  of  F igurc  9,  but  the  effect  extends  only 
a few  meters.  The  FORTRAN  computer  program 
for  the  simulation  is  given  in  the  Appendix. 

Thermal  reserve  analysis 

The  thermal  reserve  available  in  a closed  water 
body  is  often  quite  limited  and  may  prove  to  be 
the  limiting  factor  in  the  operation  of  a bubbler 
system.  In  a water  body  of  length  L,  width  B,  depth 
D,  and  average  water  temperature  T^,  the  total  heat 
content  available  for  melting  ice  is; 

(?,oui  = lbd  p,cp(r,-r„).  (28) 

As  the  bubbler  functions,  it  draws  on  this  reserve, 
lowering  the  average  water  temperature.  Since  the 


rate  of  heat  transfer  is  proportional  to 
this  decreases  correspondingly  the  heat  flux  to  the 
ice  and  the  rate  of  ice  suppression.  The  bubbler 
system  may  eventually  reach  a point  where  no  more 
ice  will  be  melted  and  ice  begins  to  re-form. 

The  computet  program  uses  the  result  of  the 
numerical  integration  of  heat  transferred  to  the  ice 
cover  and  through  open  water,  if  present,  to  cal- 
culate the  thermal  reserve  removed  from  the  water 
body  by  each  point  source  bubbler.  1 his  depletion 
of  the  ihemtd)  reserve  is  reentered  into  the  program 
by  uniformly  scaling  down  the  original  temperature 
profile  in  proportion  to  the  relative  depletion.  A 
new  impingement  temperature  is  calculated  and  the 
simulation  is  repeated.  By  appropriately  varying 
the  input  at  each  time  step,  we  may  follow  the 
evolution  of  the  suppressed  ice  covet  through  a 
weather  change  (or  season). 
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APPENDIX 


This  appendix  includes  the  FORTRAN  program  for  the  point  source  bubbler  simulation  ue- 
scribed  in  this  report.  Also  included  is  the  FORTRAN  program  for  the  line  source  bubbler  simula- 
tion. Documentation  for  the  line  source  simulation  is  presented  by  Ashton  (1974, 1978). 


POINTEiUH 

« G nSHTON  14  JUNf  1978  POINT  SOUKCf  HOBfIl  f.N  SIMULATION 

« PROGRAM  SIMULATES  OPERATION  OP  A POINT  SOURGf  PUHULPR 

« TO  MELT  ICE 

III  MI  NS  I UN  ET  A(  100)  rUUI  ( 100>  •HAT  (60) 

REAP  :t01.Hf0A.l  TA/rTUHtHWA.PELR.PEI  r 
REAP  302.AIP.AIU 
RE:AIi  303. NP 

REAP  .104.  (PAT  ( I ) . I i^l  .NP) 

301  EORMAI  ( /f  10.0) 

302  FORMAT  ( 2F 1 0 . 0 ) 

303  FORMAT  (T10) 

304  FORMAT  (FIO.O) 

4 

PRINT  401.H.UA 

401  FORMAT  (IHl.lOX.'H  * '.F6.2.'  ME  T t RH  1 OX , ' OA  = '.FIO.A.'  M3/S  ) 
PRINT  402.ETAZ.IUH 

402  FORMAT  (lOX.'ETAZ  = '.16.2.  ME  I ( RS  ' / 1 OX  . t UM  '.17.3.'  lifUC') 

PRINT  403. HUA 

403  FORMAT  (lOX.'HUA  = '.F7.2.'  U M J- PI  (i  C ) 

PRINT  404. NP 

404  FORMAT  ( lOX. 'SlMlll  AT  ION  IS  FOR  ' . I 4 . ' IiATS  ' ) 

PRINT  405 

405  FORMAT  ( 1 FIO  . 1 OX  . ' PAY  AIR  TEMP  ) 

PRINT  406. ( I .PAT ( I ) . 1=1 .Np) 

406  FORMAT  ( 1 OX . t 3 . 3X . F 7 . 2 ) 

« 

t H IS  depth  of  diffuser  (meters) 

4c  OA  IS  air  disrharae  (siiuare  meters  fer  serurMi) 

4t  FTAZ  IS  initial  ice  thickness  (meters) 

4 TWH  is  initial  water  t €*rat  ure  (doM  (^  > 

4c  HWA  IS  heat  trarisfer  coef  f i c i enf  water  to  air  <W/rtJ  lifO  > 

4 ND  IS  nijmher  of  daws  of  simulalior* 

4 DAT(I>  IS  dailw  air  temperature  (des  C ) 

4 DEI  R IS  radial  iiirreiiief*t  distance  <m) 

4 DELT  IS  time  step  (sec) 

4 AID  and  Al  U ar**  width  and  lerolth  of  the  water  hodw 

4 

FATM  * 

RHOU  * 1000. 

RHfll  * VIA, 

rr*  a A?\%, 

AM  * ?.,?A 
AKU  » 0.S4 
ALAM  * J.34F!5 
CC  * O.lf.2 
CD  = I.B3 
a * V,B07 
ANU  * 1.7VE-6 
FI  * .5.14156 

MS  TAR  --  H f F'ATM/(RhOU4G) 

HUM  * f’ATMtAl  OG<  I . M/HSTAR)/<RHnW4(:D) 

DUh  f.i)Rt(Dim) 

dUff  * (t)A440.57!))4IiUd4:».4a;4<H  + . 08 ) 4S0RT  ( R I ) 

D - <Ht().H>4r  C4(»A440. 15 
IK'H  ()WM/(.L^.4ritD4f<) 

I'RINI  407*lK.H»dUH 
IRINT  408. D 

HP  » 2. 10*AKU*(UrH»*0.55)/(  (ANII»*0.55')*(P»*0.45)  ) 

FKINT  409, HP 
NTP  06400. /PF  I T 
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mCEDlNO  PAOl  Blarac-NOT  riUOD 


40/  fWKrtA/  ( IMO^JOX.'ULH  « ' rf 8.4f 'M/S/10X» 'QWH  « 'fhB.Ar'  h3/S') 
408  fOKhOl  ilOXf'F  .FM.4f  hMKRS') 

■107  - '•Ml.J,  U,  (i:!-  DM;  C ) 

HP  BO  1-1 . I 00 
Mrt<  I ) f JO.’ 

HO  t'tlNUNlll 


DO  0t)0  U*  I.ND 

♦ tstdhlish  vjriatiofj  ot  HUl  at  DtLK* 

R - 0, 

DO  110  1*1.100 
If  (K  D) 

lOl  (JUKI)  - - 0.:*’f.*K*K/<is4D>  >*TUO 

00  ro  105 

lo;*  (tUMl)  ^ Hh*<  (D/R>i*0.4r.)*lUtl 

lot  K * KIDI  I k 

l!‘)  I'UNllNUl 

♦ Now  ralcolatf-  thic^Pnio'l  ui»rj  <iu  ltiri«  of  ire 
DO  IVT.  .f=.I.NrD 

DU  170  I *1 .100 
U liOl  M D M Ilf..  1 1'^.  1 l> 
t".  (JI  DAM  ID'  i Wt(  I ' 'Ot  I n . 'MWA) 

>.o  t U I W 
I Cf.  IM  0.0 

! I MIN  I I M4 

DMMA  ! UUM1>)  « lU  L 1 / - Miin»AL  Aei) 

M Ai  1 ) MAt  1 » ( mM  M A 

n \\  \i\\  M • 1 'I  • 1 • ! 

I .1  f f A<  M 0.0 
I UiJNIlNUl 

I Vi»  t UN  t ! Ul  l| 

l”‘.  CUNl  I4UI 

i 1 I tiU  I (-»• 

M /m  4 0 < 

DU  IV.’  I I . 1 A . 

II  <1  1A>  I I > i v.s.  I 7/..  r.'.^ 
r^.s  kl  DUI  M iMil  ♦ Ilf  Ik 
r-*  • riiNi  iNui 

Ph’lNI  410»ir»rkFDGB  ^ ^ 

410  rOkhAT  ( 1H0» lOX » ‘AF Tf K 14. DAYS  ICt  EDGf  IS  AT  R * rFA.2»  M 

4 Calculate  ihet-mal  deFletioii 

DlHERh  * FI*REDGE«R£DGf 4HWA4TWU 
NR  RFDGf/DtLK 
IF  <NR  - 100)  JOl  .I’ll. 211 
201  DO  210.  l^NR.lOt* 

Rl.R  RtDGK  ♦ Dtlk/2 

DTHFK'M  = DrUfRN  4 2.  *f  J K4ftWJ  ( I ) 

JIO  CONTINUE 
211  CONMNUF 

TWH  * 1UH4(l.  - IiTHf  kfl.MAI  D*ALW4ll*kUUW<(‘f»  ♦ 

PRINT  4II.TUM 

411  FORNAf  {iOXf  'NlU  TWH  DMi  i > 

DEtkk  * l..*DfLK 

PRINT  412.DFLRR 

412  FORhAl  (lOX.  ICt  THirKNTSS  AT  Cl  '^ND  I'M  T k HI  T[  RS  ' > 

PRINT  413.(ETACI).  I l.IOO.S' 

413  FORMAT  (lOX.lOfA.J) 

J50  CONTINUE 

END 


) 


L INF  DUD 

t G ASHTON  14  JUNE  19/fl  LINE  SOURCE  DUPDl.ER  SIMULATION 

4 

4 PROGRAM  SIMUIATES  OPERATION  OP  A LINE  SOURCE  PUPBLER 

« TO  MELT  ICF 

4 

DIMENSION  E TA(  lOOT  rOUM  1 00 )* DAT  (60) 

Rl  AD  301  .H.OA.E  TAZ.  TUH.ITNA.Df  LYfDFl  F 
Rf  AD  302 .Al  D.AI  U 
RTAD  303fND 

RIAII  304*<DAT(I  )f  1*1»ND) 


fOl 

f ORMAT 

( 7F  10.0) 

302 

FORMA r 

<;*F  10.0) 

303 

FORMAT 

( 110) 
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304  FQKHAT  (FIO.O) 

« 

PKINT  401fH»QA 

401  FOKHAF  <IHI*10X»'H  ■ HE TERS V lOX t ' OA  - 'fFlO.6#'  H?/S') 

PRINT  402fETAZ»TUH 

402  FORHAT  <10X#'ETAZ  » '»FA*2»'  HETERS  VlOXf  MWH  * '»F7.3»'  DEO  C'> 
PRINT  403»HUA 

403  FORHAT  (lOXf'HUA  *»  'fF7.2»'  U/H2-nEG  C'> 

PRINT  404»ND 

404  FORHAT  ( lOX p ' 5IHULATI0N  IS  FOR  '»14»*DAYS') 

PRINT  405 

405  FORHAT  ( IHO f 1 OX f ' DAY  AIR  TEHP') 

PRINT  40A» < 1 >DAT < I ) » 1«1 »ND> 

406  FORMAT  ( lOX » I 3- 3X » F 7 • 2 ) 

* 

t M IS  de^th  of  diffuser  (eeters) 

$ OA  IS  air  discharge  (souare  aeters  i^er  second) 

4 ETAZ  IS  initial  ice  thickness  (eeters) 

4 TUH  IS  initial  water  te»f*erature  (dee  C) 

4 NWA  IS  heat  transfer  coefficient  water  to  air  <U/H2-'D£G  C) 

4 ND  IS  nueber  of  daws  of  simulation 

4 DAT(I>  IS  dailw  air  temperature  (dee  C) 

4 DELY  IS  lateral  distance  increment  (meters) 

4 DELT  IS  lime  step  (sec) 

4 A1  B and  ALN  are  width  and  leneth  of  the  water  bodw 

4 

PATH  * 101325. 

RHOU  « 1000. 

RHOI  3 916. 

CP  » 4215. 

AKI  * 2.24 
AKU  = 0.54 
ALAH  * 3.34ES 
CC  * 0.182 
CB  « 2.14 
G * 9.807 
ANU  * 1.79E-6 
PI  » 3.14156 

HSIAR  » M ♦ PATM/(RHOU*G) 

R * ()#FO.H>4CL'4OA4*0.  15 

LilIH  = PAlMtAl  Ori<  1 . H/HfiTAR)/<RHOW4CH) 
liUM  •'.OKI  (HUH) 

(H  n ^ (UA**0.425)*HnM/(F  I4»0.254r.QRT<B)  ) 

OUH  SOkI  < 2.4F'I  )4H41ICH 
fklNT  40^«(U:M>0UM 
PkINI  4011.  H 

MR  O.  V,.«AKU4Ui:H4>0.6t  / ( AN(I440.624I<440. 38) 

FKINT  4<)'#*MH 
NTH  - 06400. /HU.  T 

107  fOKMAT  < IHO. 10X»  DCH  * ' » F 8 . 3 • ' H/S/ 1 OX • OWN  “ 'trS.S#'  H2/S'> 

108  niRMAI  llOX»  b « '»f0.3»'  HF  11  RS ' ) 

40V  MIRMAT  (IHO.'HR  « '.f«.3,’  U/H2-HEG  C.  ) 
ml  80  1 • I . 100 
M A< I > r ) TA/ 

;10  UlNlMMIt 

HO  2*-0  lli«l»N[i 

* I .t^tltsh  l.4leral  variation  of  OMI  at  DEI  Y intervals 

t ■«  0.0 

DM  MO  I Ml 00 
U <Y  It)  10U102M02 

lOl  CIUKI*  HF*4<1.|90  - 0. 1904Y4Y/(H4D)  )*TUH 
(iO  lU  103 

to.  OUKl)  > HH4(  <P/Y>440.27>4rUH 
103  Y YFHFlt 
110  MINTINIIE 

4 Ni>w  calrtilate  thicFenind  and  meltine  of  ice 

DO  195  J'lfNTD 
HO  190  1> 1*100 
IF  ( DA  r ( I D ) > 115*115*116 

115  01  ■ -IiAT(ID)/(FTA(l>/AKlkl./HUA> 
fiO  TO  117 

116  01  • 0.0 
117  iriNtlNUI 

HllltA  • (Q1  - OUKD)  4 D£LT/(RH0I4ALAH) 

F rA< I ) ■ r 1A( I ) 4 DELF  TA 


IF  (ETA<I))  1?1.122.122 

121  t IA< I ) > 0.0 

122  L'ONIINliF 
190  rONTINUF 

195  CONTINUI 

» Finrl  lev  vilav 

KFliGF  « 0.0 
DO  197  I>1.100 
IF  (F1A<I))  196>19A>197 

196  REDbF  • KEDGE  F DELY 

197  CONTINUE 
PRINT  410.ID>REDGE 

410  FORHAT  < IHOf lOXf ’AFTER  ’.lAt'DAYS  ICE  EDGE  IS  AT  Y • '»F*.2t'«') 

» Calculatv  thvraal  dvplvtion 

NR  - REDGE/DFLY 
IF  <NR  ■ tOO)  201i211>211 
201  DO  210.  I<NR>100 

RLR  > REDGE  * DELY/2 

DTHERN  • DTHERH  F 2.40UI  ( I )«DELY 

210  CONTINUE 

211  rONlINUt 

TUll  ■ TUH»(1.  - DTHCRH/(ALD*ALM*H*RHOy«CP)  > 

PRINT  411. TUFF 

411  FORHAT  (lOX.'NEU  TUH  - '.F6.3.'  DEO  C'> 

DELYY  • S.tOELY 

PRINT  412. DELYY 

412  FORHAT  (lOX.'ICE  THICKNESS  AT  CL  AND  DELYY  - ' .F5.2. 'HETERS' ) 
PRINT  413.<ETA(I>.  1*1. 100. 5) 

413  FORHAT  (10X.10F6.3) 

2S0  CONTINUE 

END 
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A facslalle  catalog  card  in  Library  of  Congreaa  MARC 
foraat  ia  reproduced  below. 
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